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ABSTRACT 

We present three-dimensional (3-D) simulations of rotationally induced line variability 
arising from complex circumstellar environment of classical T Tauri stars (CTTS) using the 
results of the 3-D magnetohydrodynamic (MHD) simulations of Romanova et al., who con- 
sidered accretion onto a CTTS with a misaligned dipole magnetic axis with respect to the 
rotational axis. The density, velocity and temperature structures of the MHD simulations are 
mapped on to the radiative transfer grid, and corresponding line source function and the ob- 
served profiles of neutral hydrogen lines (H/3, Pa/3 and Br7) are computed using the Sobolev 
escape probability method. We study the dependency of line variability on inclination angles 
(i) and magnetic axis misalignment angles (0). We find the line profiles are relatively insen- 
sitive to the details of the temperature structure of accretion funnels, but are influenced more 
by the mean temperature of the flow and its geometry. By comparing our models with the 
Pa/3 profiles of 42 CTTS observed by Folha & Emerson, we find that models with a smaller 
misaligngment angle (0 <^ 15°) are more consistent with the observations which show that 
majority of Pa/3 are rather symmetric around the line centre. For a high inclination system 
with a small dipole misalignment angle (O « 15°), only one accretion funnel (on the upper 
hemisphere) is visible to an observer at any given rotational phase. This can cause an anti- 
correlation of the line equivalent width in the blue wing (v < 0) and that in the red wing 
(v > 0) over a half of a rotational period, and a positive correlation over other half. We find 
a good overall agreement of the line variability behaviour predicted by our model and those 
from observations. 
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1 INTRODUCTION 

Classical T Tauri stars (CTTS) are thought to accrue material from 
their circumstellar discs via magnetospheric accretion (MA). In 
this paradigm, the magnetic field of the protostar truncates the 
disc at a range of radii about corotation, from where the material 
flows along the field lines and onto the photosphere. The kinetic 
powe r of the materi al is thermalized in shocks (e.g. ICamenzindl 
1990; Koenigl 1991) and is observed as a blue continuum excess 
(e.g. ICalvet & GuUbrindll998l) . while the hot m aterial within the 
funne l flows emits strongly in permitted lines (e.g. lAlencar & Basril 
2000). Observational and t heoretical aspects of the MA paradigm 
were recently reviewed bv lBouvier et alj fc007ah . 

There is now widespread observational support for the MA 
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model: CTTS are observed to have kiloga uss magnetic fields that 
are persistent over timescales of years (Qp hns-Rrull et alj 1 19991 : 
ISvmington et alJ|2005bl ; |johns-Kr ull 2007); the line profiles of hy- 
drogen are Doppler broadened to a width comparable with the 
stellar escape velocity and the line profiles often show an inverse 
P Cygni profile that arises when an accretion funnel is viewed 
against a hot spot (e.g. Edwards et alj 19941 ; lAlencar & Basril2o"o(il : 
iFolha & Emersonll200l}) ; time-dependent line profile studies indi- 
cate the the line profiles are modulated on the st ellar rotation period 
(e.g. ljohns & Basrilll995l : lBouvier et aij2007bh . 

Radiative-transfer models based on the MA paradigm are 
broadly successful in predicting line profile strengths and mor- 
phologies. Initial models were based on an simple aligned-dipole 
geometry for the accretion flow, but incorporated increasingly so- 
phistic ated physics starting from a tw o-level atom approxima- 
tion dHartmann. Hewett & Calvetlll994h . through a full statistical 
equilibrium calculation under the Sobolev approximation (SA) by 
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iMuzerolle. Calvet & Hartmanr] Jl998ah . a nd finally including SA 
plus an exact integration of the line profile dMuzerolle et aU200lh . 

The overwhelming evidence for variability, both in the 
continuum excess and the lines, of the emission from CTTS 
has lead us to examine departures from axisymmetry in the 
accretion geometry. A crude 'curtain s of accretion' model 
dSvmington. Harries & Kurosawall2005al) was able to approximate 
rather well to the observed variability, providing that the ac- 
cretion curtains had to have a relatively large azimuthal ex- 
tent. A similar, tailored model for SU Aur was also able 
to reproduce some of the observed variability characteristics 
dKurosawa. Harries & Sym ington 200^). 

The last few years has seen the publication of a series of pa- 
pers dealing with t he magnetohydrodynamical modelling of ac- 
cretion onto CTTS l|Romanova et alJl2002l ; lRomanova et al] l2003l : 
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Romanova et all [2004], hereafter R0 2. R03 and R04 respectively; 
Long. Romanova & Lovelacej|2007D . But how well do these MHD 
models agree with the observational data? The simplest test is to 
attempt to model continuum light curves from the hot spot dis- 
tributions of the models, and such simulations can reproduce the 
wide variety of observed variability. In particular for models with 
small magnetic misalignment angles the accretion funnels may ro- 
tate faster or slower than the star (meaning that the hot spots are 
not at fixed on the stellar surface), this naturally leads to the quasi- 
periodic variability that is often observed (R04). Useful though 
such comparisons are, the line profiles themselves encode much 
more detailed information on the kinematics and geometry of the 
flow. The aim of this paper is to make a quantitative comparison 
between the MHD models and spectroscopic observations by com- 
puting line profiles based on the density and velocity structure of 
the MHD calculations. 

Here we concentrate on three hydrogen transitions (H/3, 
Pa/3, and Br7). The Ha profile, although the strongest and 
most widely observed optical line, is usually contaminated by 
outflo w emission and absorption (e.g. iRei purth. Ped rosa & Lagol 
1996) and cannot be modelled by MA alone, instead requir- 
ing a hybrid code incorporating both accretion and outflow (e.g. 
Kurosawa, Harries & Symington 2006). The H/3 line is a better 
proxy for accretion, and observationally the near-IR lines show a 
high frequency of inverse P Cygni morphology, i ndicating that they 
too ar e better probes of the accretion geometry l lFolha & Emersonl 
1200 lh . 

In the following Section we describe the MHD model and the 
radiative-transfer model, and we give the results of our profile cal- 
culations in Section 3. We discuss our results in comparison with 
both earlier radiative-transfer models and observations in Section 4, 
and our conclusions are presented in Section 5. 



2 MODELS 

The basic model configuration of the central star is shown Fig.[T] 
Two important parameters in our models, the misalignment angle O 
and the inclination angle i are also shown in the figure. The former 
is defined as the angle between the rotational axis of the star (the 
z-axis) and the magnetic axis (with the dipole moment fi), while 
the latter is defined as the angle between the rotational axis and 
the direction to an observer. The rotational axes of the disc and the 
star coincide. In the following, we will describe our MHD models, 
radiative transfer models, assumed temperature structure and the 
sources of continuum radiation. 
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Figure 1. Basic model configuration. A star with its radius i?* is located at 
the origin (O) of cartesian coordinate system (x, y, z). The y-axis is into 
the paper. The rotational axis of the star coincides with z-axis. Its magnetic 
axis (with the magnetic moment fi) is inclined from the rotational axis by 
©, causing precession of magnetic axis as the star rotates. This angle will 
be referred to as the misalignment angle. The inclination angle i is defined 
as the angle, measured from z-axis, to an observer located at infinity. 




Figure 2. Example of three-dimensional simulations of MA flow for the 
magnetic axis misalignment angle = 15° . The background shows the 
iso-surface for one of density levels, red lines show sample magnetic field 
lines, and the black arrow shows the direction of the magnetic moment of 
the star. The rotational axis coincides with z-axis. 



2.1 MHD models 

Three-dimensional MHD code and model used in this paper were 
developed and described earlier in Koldoba et al. (2002), R03, and 
R04. Here we briefly discuss the main aspects of the model and 
also describe new simulations runs. 

We investigate matter flow around a rotating star with a mis- 
aligned dipole magnetic field. A star is surrounded with a dense 
cold accretion disk and a hot low-density corona above and below 
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Figure 3. Slices of density distribution (background) and sample magnetic field lines for the simulations with the magnetic axis misalignment angle = 15° 
(left), 60° (centre) and 90° (right). The X-Z plane of slices coincides with the plane defined by the rotational axis and the magnetic axis (c.f. Fig[T}. 
The red colour corresponds to the maximum of the density (p = 2 po) while the dark blue colour to the minimum density (p = 0.003 po) where po = 
4.9 X 10 — 12 g cm -3 . The units of X and Z dimensions are in R$ = 3.6 X 10 11 cm. The accretion on to the surface occurs in two streams, and the latitudinal 
location where the gas impacts on the stellar surface decreases as increases. These three density distributions (along with corresponding temperatures and 
velocities) will be used in in the subsequent radiative transfer calculations. 



it. To calculate matter flow we solve a full system of magnetohy- 
drodynamic equations (in three dimensions) using a Godunov-type 
numerical code (see Koldoba et al. 2002; R03 and R04). Equations 
are written in the coordinate system rotating with the star. A vis- 
cosity term has been added to the code with viscosity coefficient 
proportional to the a— parameter. 

The boundary conditions used here are similar to those in R03 
and R04. At the stellar surface, 'free' boundary conditions to the 
density and pressure are used. A star is treated as a perfect con- 
ductor so that the normal component of the magnetic field does 
not vary in time. There is however a 'free' condition to the az- 
imuthal component of the magnetic field: d(RB t j > )/dR = so 
that magnetic field lines have a 'freedom' to bend near the stellar 
surface. In the reference frame rotating with the star the flow ve- 
locity is adjusted such that to be parallel to the magnetic field B at 
R — R, which corresponds to a frozen-in condition. Matter falls 
to the surface of the star supersonically and most of its energy is 
expected to be radiated in the shock wave close to the surface of 
the star (e.g. lCamenzin3ll99(]| : lKoenigllll99ll : ICalvet & Gullbrind 
1998). Evolution of the radiative shock wave above the surface of 
CTTSs has been considered in detail bv lUstvugova et al.l J2006|) . In 
this paper we suggest that most of kinetic energy of the flow is con- 
verted to radiation (see R04 and Section l2"l4t . At the outer boundary 
R ~ Rmmx, free boundary conditions are taken for all variables. 

Numerical simulations have shown that the inner regions of 
the disk are disrupted by the magnetosphere of the star and matter 
flows to a star in two high-density funnels streams under some sit- 
uations. Fig.|2]shows an example of such two-stream accretion for 
a system with the misalign ment angle O = 15° (see also R04; 
Kulka rni & Romanoval 120051) . The important parameters are the 
initial densities in the disk pd and corona p c , and the initial temper- 
atures in the disk Td and corona T c . These values are determined 
at the fiducial point at the boundary between the disk and corona 
(at the inner radius of the disk). We took parameters p c — O.Olpd, 
Td = 0.01T C and obtained supersonic funnel streams since the 
sound speed is relatively low, c s ~ O.Ivk- 

The funnel flow converges towards the star and temperature 
increases due to adiabatic heating. In reality, the temperature may 
decrease due to radiative cooling. To mimic the effect of the ra- 
diative cooling we performed most of these new runs at smaller 



adiabatic index 7 = 1.1. This is the main difference of the new 
runs compared to R04 runs. In addition, we increased the magnetic 
field of the star by a factor of two compared to R04. This led to 
larger magnetospheric gaps compared to R04. 

We have adopted the following MHD input parameters from 
a typical CTTS. The mass and the radius of the star are assumed 
to be M* = 0.8 Mq and R, = 1.8 Rq respectively. The mag- 
netic field at the surface of the star (at the equator) is assumed to 
be B t = 4 x 10 3 G. The size of the simulation region corresponds 
to i? max = 40.R, = 0.34 AU. The unit of time used in the model 
is P = Pq = 1.38 d which corresponds to a period of Keplerian 
rotation at R = Ro = 3.6 x 10 11 cm which is the unit distance 
used in Fig. [3] The star rotates with period P* = 3.9 d which 
corresponds to many CTTSs (e.g. Irlerbst et al.ll2002l) . The accre- 
tion disc is stopped by the magnetosphere at the distance R t ~ 
1.67?o = 5.8 x 10 11 cm. This distance is slightly below the coro- 
tation radius of the star R COI = 2Ro = 7.2 x 10 11 cm. This situ- 
ation approximately corresponds to the rotational equilibrium state 
in w hich a star does n ot gain or lose angular momentum on average 
(see lLonget"aiT 2005). We consider the case with typical value of 
corona temperature and the disc den sity are set T c =4.5 x 10 6 K 
and po = 4.9 x 10~ 12 gem -3 (e.g. lrlartmann et aljl998l) respec- 
tively. 

R03 and R04 showed that the geometry of the accretion stream 
strongly depends on O; hence, we consider cases with a wide range 
of 6 angles, specifically 6 = 15°, 60° and 90°. Fig. [3] shows the 
density slices of the 3-D simulations on the planes defined by the 
rotation axes and the misaligned magnetic axes. The figure shows 
the accretion occurs in two streams (as mentioned earlier), and the 
flows encounter the stellar surface near the magnetic poles (see 
R04 for their relative locations). We choose the time slices of these 
MHD simulations at which the flows are semi-steady and has a sim- 
ilar accretion rate (A/ acc « 2.0 x 10 -8 Mq yr -1 ) for a comparison 
and for the radiative transfer models presented in the following sec- 
tions. This mass-accretion rate chosen here is very similar to that of 
a typical CTTS se en in the observations (e.g. Gull bringetal.1 19981: 
ICalvet et alj2004l) . 
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2.2 Radiative transfer model 



The radiative transfer code TORUS faarrieshoofll Kurosawa et al. 
2004L IKurosawa et alJ 120051 : ISvmington etai] l2005tJ : 
Kurosawa et alJ [2006) was extended to incorporate the den- 
sity, velocity and gas temperature structures from the 3-D MHD 
simulations of R04 mentioned above. The radiative transfer 
code uses the three-dimensional (3D) adaptive mesh refinement 
(AMR) grid, and it allows us an accurate mapping of the original 
MHD simulation data onto the radiative transfer grid. Although 
it is possible, we do not explore the line variability due to the 
time-dependent nature of the accretion in this paper. The aspect we 
investigate here is the variability due to the change in the viewing 
angles (of an observer) due to the rotational motion of a star and 
its magnetosphere. For this reason, we select outputs of MHD 
simulation which have (relatively) quiet stage, i.e. we choose the 
time stage of simulations which reached a (semi-) steady state. 

We emphasise that the variability associated with the time- 
dependent nature of the flow (e.g. due to instabilities) is certainly 
worth pursuing in the future since it provides us an opportunity to 
study the kinematics of the accretion flow itself; hence, it would 
provide us an additional constraint on the geometry of the magne- 
tosphere around CTTS. This is beyond the scope of this paper, but 
should be considered in a future work. 

The basic steps for computing the line variability are as fol- 
lows: (1) mapping of the MHD simulation output onto the radiative 
transfer grid, (2) the source function (S v ) calculation and (3) the 
observed flux/profile calculation as a funct ion of rotational phase . 
In the s econd step, we use the method of iKlein & CastoJ (1978) 
(see also lRvbicki & HummeJl978l : lHartmann et alJ 1994 ) in which 
the Sobolev approximation method is applied. The population of 
the bound states of hydrogen are assumed to be in statistical equi- 
librium, and the gas to be in radiative equilibrium. Our hydrogen 
atom model c onsists of 14 b ound states and a continuum. Readers 
arc referred to lHarriesI fcOOOl) for details. 

Monte Carlo radiative transfer (e.g. lHilliej[l99 ll) . under the 
Sobolev approximation, is valid when (1) a large velocity gradi- 
ent is present in the gas flow, and (2) the intrinsic line width is 
negligible compared to the Doppler broadening d ue to the bulk 
(macroscopic) motion o f gas. In our earlier models Marries 2000; 
ISvmington" e t al. 20053), this method was adopted since these 
conditions are satisfied . However, as noted and demonstrated by 
Muzerolle et al. even with a moderate mass-accretion rate 

(10~ 7 Mq yr" 1 ), Stark broadening becomes important in the op- 
tically thick Ha line. In addition, Ha is most lik ely affected by 
the wind absorption and emission components (e .g. lEdwards et al.l 
ll994l : lReipurfh et alj[l996l ; IKurosawa et alj|2006h . but the original 
MHD simulations (R03; R04) do not contain outflow/wind compo- 
nents. For these reasons we avoid modelling Ha, and concentrate 
on lines less likely to be affected by wind and Stark broadening, 
i.e. Pa/3 (mainly), Br7, and H/3 in this paper. 

When importing the MHD simulations results (Fig. [3} to the 
radiative transfer calculations, we have introduced a cut-off radius 
0~max = 7.0 x 10 11 cm = 5.5 R,) although the radial range of 
the MHD simulations extend much larger than this value. The cut- 
off radius corresponds X ~ 2.0 in Fig. [3] This indicates that the 
puffed-up density structures in the disc seen in the MHD simula- 
tions are not included in the radiative transfer models, and the prob- 
lems are rather focused on the accretion funnel parts. This is done 
to avoid the complications of adding dust opacity and finding dust 
temperature, and to avoid very high density and low temperature re- 
gions in which the source function calculation may have some dif- 



ficulty. Beyond this cut-off radius, we simply inserted a 'optically 
thick and geometrically thin disc' which is essentially a geometri- 
cal disc with no thickness and with infinitely large opacity through 
which no radiation can penetrate. This type of disc was adopted to 
imitate the obscuration of the accretion funnels and stellar surface 
by a optically thick disk. Although we understand the importance 
of including the puffed-up regions of the disc and the dust for the 
line variability problems in very high inclination system, currently 
our model are not be able handle the regions correctly. 



2.3 Temperature structure 

The temperatu r e stru cture of the magnetosphere used by 
IHartmann et al.l dl994h is computed by assuming a volumet- 
ric heating rate which is proportional to r -3 , by solving the 
energy balance of the radiative cooling rate (see Table 1 in 
IHartmann. Avrett & Edwardsll 19821) and the heating rate. However, 
in this formulation, the normalisation is arbitrary, and it has to 
be determin ed from the multiple line fitting. On the other hand, 
lMartinlfl996T) presented a self-consistent determination of the ther- 
mal structure of the in flowing gas, along the same dipole magnetic 
field geometry as in IHartmann et alJ dl982l) , by solving the heat 
equation coupled to the rate equations for hydrogen. He found that 
main heat source is adiabatic compression due to the converging 
nature of the flow, and the major contributors to the cooling pro- 
cess are bre msstrahlung radiation and line emission from Ca II and 
Mg II ions. iMuzerolle, Hartmann & Calved Jl998bh found that the 
line pr ofile models c omputed according to the temperature struc- 
ture of lMartirj (l996) do not agree w ith observations, unlike pro- 
files based on the (less self-consistent) Hart mann et alJ J 19941) tem- 
perature distribution. It is clear that the temperature structure of the 
magnetosphere is still a large source of uncertainty in the accretion 
model, and this issue should certainly be investigated more care- 
fully in the future. 

In this paper we simply consider following three cases of 
the temperature structure of the flow: (1) Hartmann-like cool- 
ing/heating case, (2) adiabatic cooling/heating only case (directly 
from the MHD calculation), and (3) isothermal case. In case of (2), 
we find that the gas temperature from the MHD is in general too 
high (no radiative cooling); hence, we introduce a scaling factor s 
which is multiplied with the original temperature (Tmhd) of the 
MHD simulations, i.e. T = s Tmhd - This is somewhat similar to 
an arbitrary normalisation constant introduced in (1). Note in all 
the models presented in this paper, s — 1.67 x 10~ 2 and the MHD 
models results with the adiabatic index 7 = 1.1 (c.f. R04). Basic 
results of the dependency of line profiles will be presented in Sec- 
tion |3.2| In the following sections, (1) and (2) will be refer to as the 
HCH and ACH temperature structures, respectively. 



2.4 The continuum sources 

We adopt stellar parameters of a typical classical T Tauri star used 
by R04 for the central continuum source to be consistent with their 
MHD simulations, i.e. its stellar radius 7?» = 1.8 Rq and its mass 
M» = 0.8 Mq (see Section 2.1). Consequently, we adopt the ef- 
fective temperature of the photosphere T p h = 4000 K and the sur- 
face gravity lo g .9* = 3.5 (cgs), and use the model atmosphere of 
iKuruc j l fl979l) as the photospheric contribution to the continuum 
flux. The parameters are summarised in TableQ] 

Additional continuum sources to be included are the hot spots 
formed by the infalling gas along the magnetic field on to the stellar 
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surface. As the gas approaches the surface, it decelerates in a strong 
shock, and is heated to ~ 10 6 K. The X-ray radiation produced in 
the shock will be a bsorbed by the gas locally, and re-emitted as op- 
tical and UV light JCalvet & Gullbrindl998l ; lGullbring et al.l2000h 
- forming the high temperature regions o n the stellar surface wit h 
which the magne t ic field inter sects. W hile Muze rolle et al.l feOOlh . 
ISvmington et alj J2005al) and I Kurosawa et aT ~J2006) used a sin- 
gle temperature model for hot spots assuming the free-falling ki- 
netic energy is thermalized in the radiating layer, and is re-emitted 
as blackbody radiation, we adopt the multi-temperature hot spot 
model of R04 in which the temperature of the hot spots is also de- 
termined by conversion of kinetic energy plus internal energy of 
infalling gas to radiation energy (as a blackbody). They consid- 
ered the position dependent matter flux crossing the inner bound- 
ary; hence, achieving a position dependent temperature of the hot 
spot, which can be written as 



Jtv| /la 
a V2 



+ w 



1/4 



(1) 



where p, p, v r and v are the pressure, the density, the radial compo- 
nent of velocity and the speed of gas/plasma, respectively. Further 
more, a is the Stefan-Boltzmann constant, and w is the specific 
enthalpy of the gas: w = 7 (p/p) (7 — 1). A typical temperature 
(area-weighted mean) of the hotspots is around 8000 K in our mod- 
els. 

We compare this temperature Th s with the effective tempera- 
ture of photosphere T p h to determine the shape of the hot spots. If 
Ths > Tph, then the location on the stellar surface is flagged as hot. 
For the hot surface, the total continuum flux is the sum of the black- 
body radiation with Th s and the flux from the model photosphere 
mentioned above. The contribution from the inflow gas is ignored 
when T hs < T ph . 



3 RESULTS 

Using Pa/3 as an example, the general characteristics of line vari- 
ability computed by the radiative transfer model will be presented 
in this section. We present models with 5 different combinations of 
the misalignment angle and inclination angle i, as summarised 
in Table [2] First, we present the continuum variability of our mod- 
els. Second, we briefly discuss the dependency of the models on 
the temperature structures (c.f. Section 12731 . Third, we discuss the 
dependency on other main input parameters (i and 0). Fourth, we 
compare line profiles for different transitions (H/3, Pa/3 and Br7). 
Finally, we will present line equivalent widths computed as a func- 
tion of rotation phase. Unless specified otherwise, we adopt the 
stellar parameters in Table Q] for a central star. The mass-accretion 
rate (M acc ) and the inner radius (Rd) of the accretion disc used in 
our models are also shown in the same table. 



3.1 Continuum Variability 

Since the line variability is closely related with the continuum vari- 
ability, we present the light curves predicated by the models be- 
fore we present the line variability results. R04 also presented the 
light curves of the 3-D MHD simulations which are basically iden- 
tical to ours. Their light curves were computed by using frequency- 
integrated flux and do not contain an optically thick disc which 
would obscure the stellar photosphere and accretion hot spots at 
high inclination angles. On the other hand, our model computes the 
light curves at a given wavelength, and contains an optically-thick 



and geometrically thin disc. Fig.|4]shows the light curves computed 
at A c = 4800 A for all the models listed in Table[2] 

Comparing the light curves from Models A, B and C which 
use the same MHD model (0 = 15° case), one can see the depen- 
dency on the inclination angle. Note that the hot spots are located 
about 30° from the poles of the rotational axis for these models 
(c.f. Figure |6j also see Fig. 7 of R04). The amplitude of the light 
curve oscillations increases as i increases up to i ~ 60° (Model B). 
Although the models with i between 0° and 60° are not shown here, 
the same trend is observed in the additional models we have run. A 
similar trend was also observed by R04 — see their Fig. 10. As i 
becomes greater than ~ 60°, the peak amplitude does not change 
greatly since a part of the second hot spot on the lower hemisphere 
becomes visible at higher inclinations. 

The visibility of the second hot spots affects the shape of the 
light curve. In Model A and B, only one hot spot (on the upper 
hemisphere) is visible during the whole or parts of rotational phase; 
hence, their light curves exhibit a single peak in one rotational 
phase. On the other hand, two hot spots are visible in Models C, 
D and E (c.f. Figures[6]and[7}. As a consequence, their light curves 
have two peaks in one rotational phase. Note that the Model C does 
not show a clear second peak, but the shape of the light curve is 
heavily affected by the presence of the second spot. 

For the models with a fixed inclination angle but with different 
misalignment angles (Models B, D and E), the oscillation ampli- 
tudes are similar to each other except for Model E that has a large 
misalignment angle (90°). The shape of the hot spots in Model E is 
extremely elongated, and has almost a belt like shape located near 
the equatorial plane (c.f. Figure|7J also see Fig. 2 of R04). The gap 
between one edge of the 'belt' like structure to the next is (~ 40°), 
causing a smaller peak-amplitude of oscillations. 

The amplitude of the light curves (Am ~ 0.4) of Models B, 
C and D are c omparable to that of A A Tau observed in V-band 
(AV w 0.5) bv lBouvier et alj J2007bl) . This suggests that the con- 
tinuum sources (Section [2.4b used in our models are quite reason- 
able. 



3.2 Dependency on temperatures structures 

We examine how the different assumed temperature structures in- 
troduced in Section 12.31 affect the line profiles. As a demonstra- 
tion, we compute Pa/3 from a system viewed at i = 60° and at the 
rotational phase t = 0.75, with a misalignment angle = 15° 
(Model B in Table|2j see also Figs. [2] and [3). 

In particular, we will compare the profiles computed with 
the HCH, ACH temperature structures (c.f. Section 12.31 ), and 
an isothermal case (Ti so = 8000 K). We choose normalisa- 
tions/scalings of the temperature structures in HCH and ACH mod- 
els such that their density-weighted mean temperatures (T mean ) of 
the gas are comparable to that of the isothermal case (~ 8000 K). 

The results are shown in Fig. [5] The profiles computed with 
the HCH and ACH are very similar to each other - they have sim- 
ilar flux in both wings, but the core flux of the HCH model is 
slightly larger than that of the ACH model. On the other hand, the 
isothermal model produces an inverse P-Cygni (IPC) profile with 
a shallower absorption wing. The difference is mainly caused by 
the warmer gas stream present in the outer part of the magneto- 
sphere (where absorption occurs) compared to that of the HCH and 
ACH models. The core flux of the isothermal case is very simi- 
lar to that of the HCH model. Overall, the profile shapes from all 
three temperature structures are very similar to each other. From 
this exercise, we find that the main physical conditions which de- 
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Table 1. Reference Model Parameters 



termines the profile shapes, for a fixed mean temperature condition, 
are the velocity field and the geometry of the funnel flows. In other 
words, as long as the mean temperature is similar, the difference in 
the temperature structure along the stream does not make a signifi- 
cant difference in line profiles at least at the temperature used here 
(8000 K). Since we find no large difference between the profiles 
from the HCH and ACH temperature structures, we adopt the latter 
(with T mcan = 8000K) in the models presented in the following 
sections unless specified otherwise. 
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Table 2. A summary of the main model parameters 



Figure 5. Comparison of Pa/3 with three different temperature structures 
of the accretion stream with © = 15° and i = 60° (Model B) at the 
rotational phase t = 0.75. Three cases considered here are: (1) the ACH 
model: the temperature structure from R04 (sol id), (2) the HCH model: tem- 
perature structure from Hartmann et al] j 19941) (dotted), and (3) an isother- 
mal (T = 8000 K) case (dashed). The normalisation of the temperature 
structure scalings are chosen such that the mass-weighted mean tempera- 
ture (T mcan ) of the flow is comparable to that of the isothermal case, i.e. 
Tmean ~ 8000 K. Overall shapes of the profiles are very similar to each 
other, but the absorption in the red wing is slightly shallower for the isother- 
mal case compared to the other two cases. The difference in the temperature 
structure in the accretion funnels does not make a significant difference in 
the profile shape as long as the mean temperature of the flows is similar. 



3.3 Dependency on inclination i 

Next we examine the dependency of the line variability on incli- 
nation angles using the MHD models with B = 15°. Fig. ^sum- 
marises the emission maps and profiles (Pa/3) for three different 
inclination angles i = 10°, 60°, and 80° (Models A, B, and C in 
Table[2]respectively), computed at four different rotational phases 
t = 0,0.25,0.5 and 1.0. The accretion onto the photosphere oc- 
curs along two streams (R04); creating two hot spots on the surface 
- one on each hemispheres. The hot spots are located about 30° 
from the rotational axis and have an 'elongated kidney-bean' shape 
(see also Fig. 2 in R04). For small inclination cases (e.g. i = 10°, 
Model A), one spot is clearly visible at all rotational phases, but for 
higher inclinations, no spot is visible at certain rotational phases. 
The visibility of spots is very important for formation of the IPC 
profile and the variability of its absorption component; hence, the 
sizes and the location distributions of spots should be understood 
for a given model. 

The good visibility of the hot spot for i — 10° (Model A) 
results in the presence of the weak absorption in the red wing of 
the model at all rotational phases; however, the weakness of the 
absorption is caused by the unfavourable alignment of spot-funnel- 
observer line of sight. For i = 60° (Model B) and 80° (Model C), 
the absorption feature in the red wing becomes most visible at the 
rotational phase when a spot is facing towards (i.e. t — 0.75) the 
observer. The largest amount of red wing absorption occurs at the 
rotational phase at which a hot spot is facing the observer and when 
the spot-funnel-observer alignment is favourable e.g. for t = 0.75 
and6 = 60°. 

The intensity maps show that most of the Pa/3 line flux contri- 
bution is from the gas in the funnels just above the hot spots (dis- 
played in purple). For a high inclination case (e.g. i = 80°), the 
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Figure 6. Pa/3 model intensity maps and the corresponding profiles computed at rotational phases t = 0.0, 0.25, 0.5 and 0.75 (from left to right) and for 
inclination angles i = 10°, 60°, and 80° (from top to bottom). The misalignment angle of the magnetic axis is fixed at © = 15° for all the models shown 
here. The intensity is shown in logarithmic scale with an arbitrary units. The physical dimension of the images are 1.4 X 10 12 cm ( ~ 11 R„) in both horizontal 
and vertical directions. The top, middle and bottom correspond to Model A, B and C in Table ^respectively. 
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Figure 7. Same as in Fig. [6] but for Models D (upper panels) and E (lower panels) (c.f. Table0. 



presence of an accretion disc greatly affects the line profiles. For 
such cases, the line of sight to the accretion funnel located below 
the equatorial plane is ob scured by the disc at all rotational phases. 
The blue-asymmetry (c.f. iFolha & Emersonl200ll) of the profile be- 
comes largest at t ~ 0.25 (for i = 80°, Model C) as little gas is 
moving toward an observer at this phase. On contrary, a half rota- 
tion later (t ~ 0.75), the profile does not become as asymmetric as 
the one at t ~ 0.25 since the self-absorption of the photons in the 
stream significantly reduces the flux in the red wing. 

The line variability behaviour of Models A, B and C are sum- 
marised in Fig. [8] The figure shows the mean spectra (phase av- 
eraged spectra), the quotient spectra (which are the original pro- 
files divided by the mean spectra) as a function of rotational phase 
(in greyscale image), and the normalized variance sp ectrum (NVS), 
which is similar to the root-mean-square spectra (c.f. Johns & Basri 
for each model. The mean spectra of the three models are 
fairly symmetric about the line centre; however, a very weak but 
noticeable amount of absorption in the red wings can be seen in the 
spectra at all i. For i — 10° and 60° cases (Models A and B), the 
flux level in their red wing becomes slightly below the continuum 
level, but the level remains above the continuum for the i = 80° 



case (Model C). Although the line equivalent width of the mean 
spectra for i = 10° is slightly smaller than that of i = 60° and 
80° cases, no significant difference is seen in three models. In ad- 
dition to the quotient spectra in Fig. [8] the phase dependent spectra 
of each model are also shown in Fig.[9]as a different representation 
of the line variability. 



The amount of the flux variability as a function of wave- 
length is summarised as the NVS. For i = 10° (Model A), sim- 
ilar levels of variations are seen in the red and blue sides. The 
NVS is double-peaked, and is symmetric around the line centre. 
This variability pattern (NVS) is very similar to that of Ha and H/3 
from t he CTTS CW Hydra (K7Ve) presented bv lAlencar & Batalhj 
d2002h. The system has a low inclination angle (i = 18° ± 10°, 
Alencar & Batalhl] l2002|) which is consistent with our Model A 
(i — 10°). The variability patterns of the red and blue sides are 
also very similar to each other as seen in the grey scale image. As i 
increases, the symmetry breaks; the fraction of peak strength on the 
red to that on the blue becomes smaller. The peak levels the NVS 
(on the red side) for i = 60° and 80° are similar (~ 10 per cent), 
but are about twice the size in the i = 10° case. 
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Figure 8. The summary of the Pa/3 spectra computed for Model A (left), Model B (centre) and Model C (right) which are viewed with inclination angles (i) 
of 10°, 60° and 80° respectively (c.f. Fig.|6). All three models have physically same accretion stream, i.e. = 15°. For each model, spectra were computed 
at 50 different rotational phases. In the bottom panels, the mean spectra of all rotational phases are shown. In the middle panels, the quotient spectra (each 
spectrum divided by the mean spectrum) are shown as greyscale images with increasing rotatio nal phases in u pward vertical direction. The greyscale image is 
scaled from 1.1 (white) to 0.9 (black). The normalized variance variance spectra NVS {johns & Basri 1995) are shown in the top panels. The mean spectra 
for i = 10° and 60° are almost symmetric about the line centre, but they do exhibit a very weak absorption in their red wings. The greyscale images shows 
that the red absorption are sometime stronger (black regions) and some time weaker (white regions) than the ones seen in the mean spectra. A little absorption 
in the red wing is seen for i = 80° case. The ratio of the amount of the variability in red wing to that in blue wing decrease as i increases. 



3.4 Dependency on misaligned angle 

Next, we compare the models with three different misalignment 
angles: 9 = 15°, 60° and 90° (Models B, D and E respectively). 
The inclination angle i is fixed at 60° for these three models. Pa/3 
profiles at rotational phase t — 0, 0.25, 0.5 and 0.75 along with 
the corresponding spatial intensity maps are shown in Figs. [6] (for 
Model B) and [7] (for Models D and E). Although the shapes of 
the funnel flows are different, the accretion still occurs in two 
streams for all these models, causing two hot spots on the stellar 
surface. The width of the stream becomes wider as 6 increases, and 
consequently the azimuthal extent of the hot spots also becomes 
wider. The latitudinal position of the hot spots becomes lower as 
the misalignment angle increases (c.f. Fig. [3}. For = 90° case 
(Model E), the two wide and thin funnels are located almost on 
equatorial plane (c.f. Fig. [3}, and so are the hot spots. See R04 for 
larger and clearer depiction of the hot spot geometries and their 
physical properties. Interestingly, similar equatorial flows and the 
shape of hot spots are found in a quadrupole magnetic field accre- 
tion model of Lo ng et al. U2007t) with a small misalignment angle. 

The IPC profiles are found at t = 0.75 for 6 = 15° (Model B) 
and 9 = 60° (Model D), but not for 6 = 90° (Model E). 
The nearly equatorial accretion flows in Model D do not have a 
favourable spot-stream-observer line of sight which is essential for 
the formation for of the IPC profile. The double-peaked profiles are 
seen for larger models (Models D and E). The splitting of the 
peaks are caused by the rotational motion of the magnetosphere. 

Figure [10] shows the summary of the line variability (Pa/3) 
for Models B, D and E. The mean line profile becomes wider as 
increases. The separation of two peaks in the mean profile is 
largest for = 90° . The peak levels of the NVS are similar for 
all three models, but the location of the peak(s) is different. While 
the amount of variability is largest at the line centre for Model E 



(0 = 90°), it is largest in the blue wing (v ~ 40 km s -1 ) for Mod- 
els B and D. We note that variability pattern seen in M odel E resem- 
bles tha t of H/3 from the CTTS AA Tau observed bv lBouvier et al.l 
(2007b) although their inclination angle (i w 75°) is different from 
the model presented here (i — 60°). See also Fig. [9] for the sum- 
mary of the phase dependent spectra. 



3.5 Comparison of different lines 

Once again using the models with = 15° (Models A, B and 
C) as examples, we demonstrate the difference in the line shapes 
among different hydrogen lines: Pa/3, Br7 and H/3. The summary 
of the comparison spectra is shown in Figure QT] Although they 
differ in strength, the overall dependency of the line profile shape 
on the inclination angle is very similar in these three transitions. 
The relative line strength slightly decreases from H/3 to Pa/3 then 
to Br7. The line cores are narrower for higher inclination angle 
models. The strength of the blue wings are similar for i = 10° and 
60°, but much weaker in the higher inclination model i = 80° since 
most of the flow moving toward an observer (blueshifted flows) is 
eclipsed by the stellar surface at a higher inclination (c.f. Fig.[6jl. 

The line variability behaviour for each line computed at i = 
60° is summarised in Figure[T2] All three lines show a very similar 
variability pattern in both the NVS spectra and the grey scale quo- 
tient spectral images. Interestingly, the peak level of NVS for H/3 
is about 4 times larger than those of Pa/3 and Br7. The difference 
in the size of variablity can be explained by a much larger spatial 
extent of the line emission regions for H/3 compared to Pa/3 and 
Br7. 
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Figure 11. Comparison of H/3, Pa/3 and Bry for i = 10° (Model A: solid), 
60° (Model B: dotted) and 80° (Model C: dashed) from left to right re- 
spectively. All the profiles are computed at the rotational phase of 0.75 (c.f. 
Fig. |6j. Overall dependency on the inclination angles are similar for all 
three transitions although the relative strength slightly decreases from H/3 
to Pa/3 then to Bi'7. The line cores are narrower for higher inclination angel 
models. The strength of the blue wings are similar for i = 10° and 60°, 
but much weaker in the higher inclination model i = 80° since the most 
of the flow moving toward an observer (blueshifted flows) is eclipsed by 
the stellar surface in a higher inclination (c.f. Fig.|6j. The depth of the red 
wing absorption is largest for i = 60° cases where the aligment of a hot 
spot-infall stream- observer is near optimum. 



Figure 9. The time-series spectra of Pa/3 from Models A through E (from 
left to right,) are shown for 10 different rotational phase. From the top to 
the bottom the phase changes from to 1. Each profile is separated by the 
rotational phase of ~ 0.1, and shifted upward by 1.0 as the rotational phase 
increases for clarity. 



3.6 Line equivalent width 

The phase dependency of the line EW (Pa/3) from each model in 
Table|2]is shown in Figure[T3] We compute three different types of 
EWs: (1) by using all velocity bins in the model profile (from —500 
to 500 km s _1 ), (2) by using only negative velocity bins (the blue 
wing: -500 to km s _1 ), and (3) by using only positive velocity 
bins (the red wing: to 500 km s^ 1 ). Simil ar methods are often 
used in time-series sp ectra observations (e.g. I Johns & Basrill 19951 ; 
I Kurosawa et alJr2005h to examine the gas kinematics of the flow. 
Note that here we use the convention of the sign for the line EWs 
as negative when the flux is below the continuum and positive for 
vice versa. 

For the model with G = 15° and a mid to low inclination an- 
gle (i.e. Models A and B), the EW curves show only one minima 
(at phase t — 0.75) in one rotational period. On the other hand, 
with the same physical model (0 = 15°) but with a high inclina- 
tion angle (i.e. i — 80° as in Model C), the shapes of the EWs 
are affected by the presence of a second local minima at t — 0.25, 
especially in the red wing (positive velocity bins) EW curve which 
clearly shows two local minima in one rotational period. The lo- 
cal minima at t = 0.75 is caused by the maximum continuum flux 
(c.f. Fig. [4} contribution from the hot spot in the upper hemisphere 
(c.f. Fig.[6]l. The second local minima at t = 0.25 is caused by the 
combination of the following two reasons. First, at the high inclina- 
tion, only one accretion arm is visible to an observer. Furthermore, 
at this rotational phase, the flow component which is moving away 
from the observer almost completely disappears because of the ori- 
entation of the remaining upper accretion arm (c.f. Fig. [6]l hence 
causing the minimum flux in the red wing, which corresponds to 
the minimum EW of the red wing. Secondly, the visibility of the 



second hot spot (c.f. Section [3~Tt from the lower hemisphere at this 
rotational phase causes a slight increase in the continuum flux as 
also seen in the light curve in Fig. [4] This slightly weakens the line 
strength even further. 

For larger misalignment models (Models D and E), we ob- 
serve two local minima in one rotational period. As mentioned ear- 
lier, the two local minima are seen because of the visibility of two 
hot spots in these models at this inclination (i — 60°). Remember 
that the continuum light curves of these models (Fig. [4} also show 
two local minima in one rotational phase. We note that the position 
of the local minima in the EW for these two models coincide well 
with the position of the maxima (or local maxima) of the contin- 
uum light curves in Fig. [4] This clearly demonstrates that the hot 
spot visibility greatly influences the line EW variability. To quan- 
tify the statement above , we have compu ted the Pearson's correla- 
tion coefficients (r, c.f. lBevington|[l969l) for the continuum light 
curves in Fig. [4] and the EW curves in Fig. [T3] for each model. The 
results are r = 0.67, 0.92, 0.95, 0.88, 0.99 for Models A, B, C, D 
and E respectively — indicating fairly good correlations between 
the continuum and the EW variations. 

The EW of the red wing and that of the blue wing of each 
model are well correlated with each other except for Model C 
which has a very high inclination (i — 80°). The figure shows 
that during half of the rotational period (between t = and 0.5), 
the EW curves are anti-correlated each other. This anti-correlation 
is naturally caused by the fact that only one accretion arm is visi- 
ble to an observer at this high inclination. When the arm is on the 
opposite side of the star (when the star is in between the observer 
and the arm), the EW of the blue wing should be maximum and 
that of the red wing should be minimum. When the arm is on near 
side of the observer, the opposite should occur if the hot spots are 
not visible to the observer. As explained earlier, the high visibility 
of the hot spot on the upper hemisphere causes the increase in the 
continuum flux resulting in lower EW in the red wing for the phase 
between t — 0.5 and 1.0, causing the correlation between the EW 
of the red wing and the blue wing (instead of an anti-correlation). 

In the time-series spectra observation of the CTTS SU Aur, 
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Figure 10. Same as in Fig. [8] but for different magnetic misalignment angles: = 15° (left, Model B), 60° (centre, Model D), and 90° (right, Model E). 
The inclination is fixed at i = 60° for these models. The mean spectrum (bottom panels) becomes wider and more asymmetric (around the line centres) as 
increases. A weak absorption in the red wing are seen in the mean spectra for = 15° and 60° cases, but no clear red absorption is seen in = 90° case. 
The variability patterns (middle and top panels) for three cases are quite distinct from one another. 



John s & Basril Jl995h found an anti-correlation between the blue 
wing EW of Ha and the red wing EW of H/3. They discussed that 
the anti-correlation might be explained by a tilted dipole magneto- 
sphere with outflow components from the disc-magnetosphere in- 
teraction region. Although there is no wind/outflow component in 
our model, the mechanism producing the anti-correlation seen in 
Model C may be also an additional reason for the observed anti- 
correlation in SU Aur. We also note that the high inclination which 
is required for this mechanis m is consistent with the high incli- 
nation angle of SU A ur (e.g. iMuzerolle et alj|2003l ; IUnruhll2004 
iKurosawa et"aHl2005h. In time-series spect roscopic observation of 
Pa/? from SU Aur, IKurosawa et alj b005l) also found the EW of 
the red wing anti-correlates with that of the blue wing at some ro- 
tational phases, but they correlate in other phases. Once again this 
behaviour is quite similar to that one seen our Model C. 



4 DISCUSSION 

4.1 Comparison with observations 

Although not based on simultaneous observations, 
iFolha & EmersorJ fcOOll) presented 42 Pa/3 and 30 Bi'7 pro- 
files of CTTS and weak T Tauri stars (WTTS) mainly in the 
Taurus-Auriga complex. They found that 53 per cent of Pa/3 pro- 
files have shapes sym metric around the line centre (Type I profile: 
iReipurth et ail [l 996) and 34 per cent of them have IPC profiles. 
Similarly for Br7, they found 72 and 20 per cent for Type I and IPC 
profiles respectively. They also found no blueshifted absorption 
components in both Pa/3 and Br/-y except in one object (CW Tau) 
and only in Pa/3. 

All of our models in Figs.[6]and|7](except for Model E) show 
the IPC at some rotational phases. For Br7, in general we found 
that the IPC feature (the redshifted absorption component) is much 
weaker than that seen in the Pa/3. For example, Fig. [TT] shows the 
flux level in the red wing remains above the continuum for mid to 
low inclination systems, and it is below the continuum level only 



for a very high inclination system. Using the time-series model 
spectra (which are equally sampled in rotational phase) shown in 
Figs. [8] [TO] and [J_2] we have computed the fraction of the occur- 
rence of IPC profiles in each model. The results are summarised 
in Table [3] along with the fractio n of the IPC profiles in t he Pa/3 
and Bpy from the observations of IFolha & Emersonl |2001). For a 
fixed O value (Models A, B and C), the fraction of the IPC profiles 
decreases as the inclination angle i increases for Pa/3. The largest 
inclination model (i — 80°, Model C) shows the fraction of IPC 
profiles is 37 per cent which is closest to the observed value of 
34 per cent. On the other hand, no such clear trend is seen for 
the Pa/3 models with the fixed i = 60° value but for different 

values (i.e. Models B, D and E). At this mid-inclination angle 

1 = 60° (cos i = 0.5), the fraction of IPC profiles (for Pa/3) is 50, 
50 and per cent for Models B, D and E respectively. The frac- 
tions for Models B and D are much larger than that of the observa- 
tions, and that for Model E is much smaller. To compare the models 
more strictly with the IPC fractions of the observations, one needs 
enough sampling of the model profiles in both rotational phase and 
inclination angle. Although we have enough phase sampling points 
(50 angles), we lack a full sampling of inclination angle (3 angles); 
hence, the direct comparison is difficult here. 

The full width half maxima (FWHM) of the mean (phase- 
averaged) Pa/3 and Br7 spectra of our models (Figs. [8l [Tol and [T2t 
are also summarised in Tablef3]along with the mean FW HM of the 
observed Pa/3 and Bpy profiles from Folha & Emerson (200l]). The 
table clearly shows that in all models presented here, the predicted 
FWHM underestimate the observed values. The discrepancy may 
be caused by one or a combination of the following: (1) the rela- 
tively small size of the magnetosphere in the MHD models in which 
the gas velocities just before reaching the stellar surface is rela- 
tively small ~ 200 km s~ x , (2) the spatial widths of the accretion 
funnels are underestimated, (3) other line b roadening mechanisms 
(besides the Doppler) may be important (e.g. IMuzerolle et alj200lt 
IKurosawa et alj20 06). and (4) the temperature of funnel streams is 
underestimated. An implication for possible cause (2) can be seen 
in the tendency of increasing FWHM values as increasing O val- 
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Figure 12. Same as in Fig.[8] but for Models B (i = 60° and = 15°) with different lines: H/3 (left), Pa/3 (centre) and Br7 (right). The line strength of the 
mean spectra (bottom panels) decreases from H/3 to Pa/3, and then to Br7. All lines show a similar variability pattern as seen in their quotient spectra greyscale 
images (middle panels); however, they differ in their magnitudes (top panels). 



ues. Note that the widths of the accretion streams become wider 
for a larger O model. As mentioned before (Section [2.3t . the nor- 
malisation of the temperature structures used in our models is quite 
arbitrary, and the temperatures used here may underestimate that of 
a typical CTTS. 

We also find that among our models only O — 15° models 
(Models A, B and C) show rather symmetric Type I profiles. The 
larger O models (Models D and E) show profiles split near the line 
centres, and have an appearance of a double-peaked profile. The 
difference is mainly caused by the geometry of the accretion fun- 
nels and the locations of the hot spots where the accretion funnels 
meet the stellar surface. In Models D and E with O = 60° and 90° 
respectively, the accretion funnels are much wider and the hot spot 
latitudes are much lower compared to those of the O = 15° mod- 
els. This would allow an observer to see near the base of the accre- 
tion funnels from both hemisphere for the larger misalignment an- 
gle models (c.f. Fig.[7j>. This and the rotational motion of the mag- 
netosphere causes a double-peaked profiles. On the other hand, in 
the low misalignment angle models, the base of only one accretion 
funnel is visible to an observer (c.f. Fig. [6jl. Based on the compar- 
ison of the line profile shapes, we find that a model with a smaller 
misalignm ent angle (6 ~ 15°) is more consistent with the obser- 
vations of Folha & Emerson] d200lf). Note that th e double-peaked 
profile shapes (Type II-R, c.f. Reipurth et al.lll996[) seen in Mod- 
els D and E are not en tirely inconsistent with the observations, as 
Folha & E merson! d200 lb found such Pa/3 profiles in a few systems. 

The hot spot surface coverage fractions in our models (in 
Figs.[6]and[7]> are approximately 7 per cent (for Models A, B and 
C), 9 per cent (for Model D) and 6 per cent (for Model E). These 
values seem rather large compared to the hot spot coverage in- 
dicated by observations. For example, ICalve t & Gullbring (1998) 
found that for the majority of CTTS, the filling factor is / = 
0.001 — 0.01, and the corresponding surface cover age of the spots 
is onl y 0. 1-1 per cent of the stellar surface (see also Gullbri ng et al.l 
2000). Indeed, these size s are i n agreement with those derived 



Model/Obs. 


FWHM 

(km s _1 ) 


IPC fraction 
(per cent) 


Observation (Pa/3) 


204 ± 12 


34 


Model A (Pa/3) 


42 


100 


Model B (Pa/3) 


70 


50 


Model C (Pa/3) 


70 


37 


Model D (Pa/3) 


113 


50 


Model E (Pa/3) 


129 





Observation (Bry) 


207 ± 26 


20 


Model B (Bry) 


70 


31 



by Valenti. Basri. & Johns! 0_993) using s pectrophotometr i c obser- 
vations of 'blue continuum.' However, Muzer olle et al. I dl998bl) 



Table 3. Comparison with the observation of Folha & Emerson (2001). See 
Table|2]for model parameters. 



pointed out that the filling factor / seems to be much larger at 
longer wavelengths, and it was suggested that possibly the filling 
factor depends on wavelength. Here we note that the hot spots ob- 
tained in the 3-D MHD simulations are strongly inhomogeneous 
i.e. the kinetic energy per unit area (matter flux) is much larger in 
the centre of the spot c ompared to peripheral regions of the spot 
(c.f. see Figs. 2 and 3 in lRomanova et al]|2004l) . 

In addition, the line strength and shape of the models are not 
expected to be very sensitive to a hot spot coverage, provided that 
the mass-accretion rate and hence the hot spot luminosity is un- 
affected by the change in the size of hot spot itself. However, if 
the hot spot coverage becomes much larger (e.g. 30 per cent), then 
the geometry of the magnetosphere itself must change significantly. 
This would result in a significant change in the line shapes. In gen- 
eral, the lines are more sensitive to mass-accretion rates and con- 
sequently to hot spot luminosities. A possible effect of a reduced 
hot spot coverage area in the model would be a smaller chance 
of producing the IPC profiles. In turn, this would give us a bet- 
ter agreement with the IPC fraction seen in the observations of 
Folh a & Emerson! ( feOOll) . as discussed earlier (c.f. Table[3}- 
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0.5 1 1.5 2 
Rotational Phase 

Figure 13. Equivalent width (EW) variations as a function of rotational 
phase for the six models (A through F) in Table [2] For each model, the to- 
tal EW (solid), the EW in the red wing (dotted), and that in the blue wing 
(dashed) are separately computed. Except for Model C, the temporal varia- 
tions of the EW in the red wing and that in the blue wing are well correlated. 
For Model C, they are anti-correlated between phase t = and 0.5 (see text 
for explanation). 

4.2 Intrinsic variability 

In the previous section (Section [3}, we presented the 'rotationally 
modulated' variability of the spectra based on the MHD simula- 
tions. We observed that the variability associated with the absorp- 
tion in the red wing reaches up to 30-40 per cent (c.f. Figs. [8l ITOl 
and!12t. Here we briefly examine whether such rotationally mod- 
ulated variability can be observed when, in reality, the accretion 
rate through the funnels to the stellar surface may also vary in time. 
The MHD simulations have shown that in many cases the mass- 
accretion rate (M) to the stellar surface reaches quasi-stationary 
states, but the mass flux continues to vary on a time-scale of 3-5 d, 
which is similar to the rotationally modulated time-scale (~ 4 d). 
Thus the line profile variability induced by the stellar rotation may 
be complicated by the intrinsic variation of the accretion rate. How- 
ever, the MHD simulations show that in one of the main cases 
(O = 15° and 7 = 1.1), the accretion rate M is almost constant, 
and varies only at the level of 5-6 per cent. In this case, we ex- 
pect that the effect of the intrinsic variability on the line variability 
is relatively small, and the results presented in the previous sec- 
tion are not affected significantly. On the other hand, in case for 



9 = 15° and 7 = 5/3 (c.f. lRomanova et alj|2003l,l2004l) . we find 
that the variation of M occurs at the level of 10-30 per cent. We 
expect in this case, the line profile variability to be affected by the 
change in M, which occurs on a similar time scale as the rotational 
period but varies rather stochastically. The line variability will be a 
mixture of periodic signals from stellar rotation and rather random 
signals from the change in M. However, we expect that in this case 
long simulation runs will help to separate the contribution from this 
large but random component (caused by the change in M) from the 
total line variability. 



4.3 The Sobolev approximation 

Although not explicitly shown here, the funnel flow velocities near 
the inner edge of the accretion disc are rather slow and subsonic. 
For example in Models A, B, and C (© = 15° cases), the flow 
ve locities are very s imilar to the one shown in Fig. 5 (right panel) 
of Romano va et all d2004l) in which the Mach number of the gas 
flow in the middle of the funnel flow is plotted as a function of dis- 
tance from the star (along the stream line). The plot shows that most 
of the flow is in fact supersonic except for the gas near the disc. 
The Sobolev Approximation used in the radiative transfer models 
is not quite valid in the subsonic part of the flow; however, most of 
the line emission occurs in the funnel flow near the stellar surface 
(c.f. Figs. [6] and [7}, and the line of sight to the emission region 
does not pass through the slow moving part of the funnel flow near 
the inner edge of accretion disc, except for a very high inclination 
case. For this reason, the Sobolev approximation used in the line 
profile models presented here should be a reasonable assumption, 
provided that the intrinsic line width is negligible compared to the 
Doppler broadening due to the bulk (macroscopic) motion of gas. 

In Section l4~Tl we noted that the possibility of additional line 
broadening mechanisms besides the the Doppler may be impor- 
tant for explaining much larger FWH Ms seen in the ob s erved Pa/3, 
compared to th o se fro m our models. iMuzerolle et al. I J200lh and 
iKurosawa et alj l[2006j) considered Stark broadening which is im- 
portant in optically thick lines e.g. Ha, but less likely important 
for Pa/3 and Br7 lines. Another possible cause of line broadening 
is due to large turbulent motion of plasma which c ould be caused 
by Al fven waves excited in the magnetosphere (c.f. I Johns & Basril 
Il995l) . If the magnitude of turbulent velocity is comparable to that 
of the local thermal velocity, the assumption of the Sobolev ap- 
proximation becomes invalid. The MHD simulations presented in 
this paper do not show such turbulent motions; hence, we have sim- 
ply adopted the Sobolev approximation in the work presented here. 
On the other hand, if large turbulence is found in the MHD simu- 
lations, the radiative transfer models should abandon the Sobolev 
approximation and us ed an alternative metho d e. g. the ray-by-ray 
integra tion method in IMuzerolle et al] d200ll) and lKurosawa et all 
J2006h . 



4.4 Comparison with earlier models 

Although very successful in e xplaining many CTTS emission 
line features, the MA model of IMuzerolle et al] d200ll) are axi- 
symmetric (2.5-D); hence, they cannot predict the line variability 
associated with stellar rotation. On the other hand, their model still 
can predict line variability c aused by a non-constant mass-accretion 
rate. Kurosawa et al. (2006) introd uced the disc wind - MA hybrid 
mo del by combining the m odels of K nigge. Woods & Drew! ( fl995T) 
an d Hart mannetaiT ( ll994h i.e. including both outflows and inflows 
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from and to the CTTS. The model was successful in reproducing 
the wide variety of Ha profile shapes seen in observations, and it 
demonstrated the importance of a wind or outflow component in 
determining the profiles shapes of Ha. Their models were also axi- 

symmetric, and did not consider line varia bility. 

The same dipo l e MA model as in Hart mann et alj (l994) 
and iMuzerolle et al.1 fcOOll) are used in the 3-D M onte Carlo 
radiative transfer mo del of ISvm ington et al.l (2005a) (see also 
iKurosawa et alj |2005). They have modified the magnetosphere by 
removing the flow within some azimuthal angle ranges, and then 
displaced the magnetic poles from the rotational poles to imitate 
the flow geometry found in the MHD simulation of R03 and R04, to 
simulate a possible line variability for i — 60° case. The misalign- 
ment angle used in this model (their A30 model) is rather small 
(10°), and it is very similar to our = 15° models. Although the 
line strength is much weaker than our Model B (which also uses 
i — 60°), the variability pattern seen in their RMS spectra and the 
quotient spectra image (see their Figure 9) are very similar to ours 
(Fig. [8). The difference in the line strength is caused by the dif- 
ference in the adopted mass-accretion rate and the mass-weighted 
mean temperature between the two models. 



5 CONCLUSIONS 

We have presented a series of 3-D Monte Carlo radiative trans- 
fer calculations of the line variability from CTTS induced by 
stellar rotation using the resu lts of the 3-D MHD simulations of 
iRomanova etal] (2003, 2004) who considered the accretion onto 
the CTTS with a misaligned dipole axis with respect to the rota- 
tional axis. Our model here does not include an outflow component, 
and is restricted to the case of accretion flow through the magne- 
tosphere. Main objectives here were firstly to examine whether the 
3-D MHD simulations would be able to reproduce the types of line 
profiles and the line variability seen in observations, and secondly 
to examine how the line variability depends on basic physical pa- 
rameters. In the following, we will summarise our main findings 
though this investigation. 

(1) The detail of the temperature structure along the accretion 
funnel flow does not affect the line shapes greatly (Section 13.21 ). 
For a fixed (density- weighted) mean temperature (Tinean) of the ac- 
cretion flow, we have considered three different temperature struc- 
tures (HCH, ACH and isothermal c.f. Section l2~3t . and found the 
resulting Pa/3 profiles are very similar to each other, at least with 
Tmcan = 8000 K. The line shapes are more sensitive to the incli- 
nation of the systems, the geometry and hence the velocity field of 
the accretion flow. 

(2) By comparing our models with the a tlas of observed 
Pa/9 and Br7 profiles of iFolha & Emersor] |2001), we found mod- 
els with a smaller misalignment angle (e.g. Q = 15°) produce 
rather symmetric profiles which are seen in the majority of the ob- 
served sample. The models with mid to high misalignment angles 
produces double-peaked profiles (Type II-R: Reipurth et al. 1996) 
which are very rare in the sample of lFolha & Emerson ( 2001 ). This 
may also suggest that the majority of the CTTS have a rather small 
misalignment angle. 

(3) We find that the line equivalent width variability is closely 
related with the visibility of the hot spots on the stellar surface 
(Sectio rl3~6l . For a high inclination system with a small dipole mis- 
alignment angle (e.g. Model C), only one accretion of funnel (on 
the upper hemisphere) is visible to an observer at any given rota- 
tional phase; hence, the anti-correlation of the EW in the blue wing 



and that in the red wing can occur during half of a rotational pe- 
riod, specifically when the hot spot is not visible to the observer 

(Fig.[H). 

(4) Ba sed on our line profile models , we find that the MHD 
models of Ro manova et al. f J2003I . 120041) are capable of repro- 
ducing line profil e variability behaviour similar to those seen in 
observations (e.g. Alencar & Batalha|[2002l ; lKurosawa et al1l2005l : 
iBouvier et ai1 l2007b) although the original temperature predicted 
by the MHD models had to be lowered by an arbitrary scaling factor 
in the radiative transfer calculation due to lack of a proper cooling 
mechanism in the MHD model (c.f. Sections [2.3| l. 

Despite of the relatively good agreement between our mod- 
els and observed line variability, several issues still remain. First, 
the assumption of the Sobolev approximation used in the radiative 
transfer calculation may not be valid in the funnel flow near the 
inner edge of the accretion disc where the speed of flow is substan- 
tially subsonic. Relaxing this assumption may be more important 
in the observed flux calculatio n than in the source f unction calcula- 
tion according to the finding o flBastian et al. I [Hi) who compared 
line profiles computed in 1-D (spherical geometry), using three dif- 
ferent radiative transfer methods including the Sobolev approxima- 
tion method. Detailed investigation of the validity of the Sobolev 
approximation in the context of the magnetospheric accretion is 
underway, and shall be presented in a future paper (Harries et al., 
in preparation). 

Second, the model presented here does not include a wind 
component. Although the effect of the wind on the line profile 
shapes are expected to be small in Pa/3 and Br7, it could poten- 
tially have significant effect s on H/3 and Ha (e.g., Reipu rth" et al.l 
ll996l ; lAlencar & Basri 2000). Strong wind absorpti on components 
are al s o seen the optically thick He I A 10830 lin e tedwards et al.l 
120031 : iDupree et alj 120051 ; lEdwards et al] 120061) . These spectro- 
scopi c observations combined with models (e.g ., iMa tt & Pudritz 



2005, 20071; iFerreira. Dougados & Cabritl 120061 : IKurosawa etal 



2006; Kwan. Edwards & Fi scher 2007) provide us with opportuni- 
ties to constrain little known physical properties (e.g. geometry and 
temperature) of the wind/outflow in sub-AU scales. However, as the 
MHD models used in this paper do not include the outflow, we are 
not be able to explore the possible effect of the wind contribution 
to line variability. In principle, our radiative transfer model can be 
applied to a system which includes both the MA flow and the wind 
once a result of such MHD calculations becomes available. 

Strict tests of our models should be performed by quantita- 
tive fitting of the time-series observations of m ultiple lines (e.g., 
Ilohns & Basrlll 19951 ; IUmir]|2004lBo"uvier et al.ll2007bh with good 
rotational phase coverages. The different lines are formed in dif- 
ferent volumes within the funnels and give us additional informa- 
tion on the flow kinematics. The spectra from different rotational 
phases provides information on how the magnetosphere changes in 
azimuthal direction. Once the validity of the model is established, 
it can be used to constrain the complex geometry of the MA flows 
around CTTS. It would be very interesting to compare the geom- 
etry constrained by this method to the ones determined from the 
field extrapolation method (e.g. ljardine. Collier Cameron & Donatil 
120021 ; iGregorv et al] [2006 ) using the stellar surface magnetic field 
information via the Zeeman-Doppler imaging technique, obtained 
with a modern stellar spectropolari meter such as ESPaDOnS 
jPonati et alJl2006l ; lDonati et al.ll2007h . 

While our underlying assumption of the magnetic field ge- 
ometry is pure di pole, a more complex field geometry, such as in 
Long et al. (2007) who considered the combination of a dipole and 
a quadrupole fields in their MHD models, should be also considered 
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in the line variability calculations. For example , some CTTS are 
known to posses non-dipole magnetic fields (c.f . Johns-Krull et al.l 



1 19991 : IValenti & Johns-Krullll20ol iDonati et alJl2007h . Finally, we 
shall also consider the case with a non-constant mass-accretion 
rate model. For example, a rec ent MHD calculation b y Kulkarni 
& Romanova in preparation (see Romanova et al. 2007) has shown 
that variable mass-accretion rates could occur due to instabilities 
(Rayleigh-Taylor and Kelvin-Helmholtz) in the MA flow. These 
add another level of complexity in the interpretation of observed 
line profiles, and their variability. 
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